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ABSTRACT 

We extend the method introduced by Cinzano et al. (2000a) to map the artificial 
sky brightness in large territories from DMSP satellite data, in order to map the naked 
eye star visibility and telescopic limiting magnitudes. For these purposes we take into 
account the altitude of each land area from GTOPO30 world elevation data, the 
natural sky brightness in the chosen sky direction, based on Garstang modelling, the 
eye capability with naked eye or a telescope, based on the Schaefer (1990) and Garstang 
(2000b) approach, and the stellar extinction in the visual photometric band. For near 
zenith sky directions we also take into account screening by terrain elevation. Maps of 
naked eye star visibility and telescopic limiting magnitudes are useful to quantify the 
capability of the population to perceive our Universe, to evaluate the future evolution, 
to make cross correlations with statistical parameters and to recognize areas where 
astronomical observations or popularisation can still acceptably be made. We present, 
as an application, maps of naked eye star visibility and total sky brightness in V band 
in Europe at the zenith with a resolution of approximately 1 km. 

Key words: atmospheric effects - site testing - scattering - light pollution 



1 INTRODUCTION 

The recent availability of high spatial resolution radiance 
calibrated night-time satellite images of the Earth (Elvidge 
et al. 1999) allows one to obtain quantitative information 
on the upward light flux emitted from almost all countries 
around the World (e.g. Isobe & Hamamura 1998) bypassing 
problems arising when using population data to estimate 
light pollution: (i) census data are not available everywhere, 
(ii) they are not updated frequently, (iii) they are based on 
city lists and do not provide spatially explicit detail of the 
population geographical distribution, (iv) they do not well 
represent some polluting sources, like e.g. industrial areas, 
harbours and airports, (v) the upward emission per capita 
of a given city can deviate from the average and geographic 
gradients could exist. 

In the last 16 years Roy Garstang has been carrying on 
a strong modelling effort (Garstang 1984, 1986, 1987, 1988, 
1989a, 1989b, 1989c, 1991a, 1991b, 1991c, 1992, 1993, 2000a) 
to develop an accurate technique to evaluate the night sky 
brightness produced by the upward light flux. It avoids re- 
sorting to empirical or semi-empirical formulae which do not 
allow detailed relations with the atmospheric conditions, 
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choice of the direction of observation and accounting for 
Earth curvature, even if they are of invaluable utility for sim- 
ple estimates (e.g. Walker 1973; Treanor 1973; Berry 1976; 
Garstang 1991b). 

Cinzano et al. (2000a, hereafter Paper I) applied 
Garstang models to DMSP satellite data to produce detailed 
maps of the artificial night sky brightness across large ter- 
ritories opening the way to a quantitative analysis at global 
scale of the entire Earth (Cinzano et al. in prep.) and, joined 
to an even more continuous observation of the Earth made 
by DMSP satellites, to the prediction of the future evolution 
(Cinzano et al. 2000b; Cinzano et al. in prep.). 

Both a comprehensive study of the effects of the increase 
of light levels in the night environment over their natural 
condition produced by wasted light and the evaluation of 
the effectiveness of laws, standard rules and ordinances to 
protect the environment and the capability of mankind to 
perceive the universe, require more than maps of artificial 
sky brightness at sea level. Such maps, being free from ele- 
vation's effects, are useful for a detailed knowledge and com- 
parison of the pollution levels across large territories and the 
recognition of most polluted areas or more polluting cities. 
They are also useful for the identification of dark areas and 
potential observatory sites. However they allow only in an 
approximate way a quantitative evaluation of the capabil- 
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ity of the population to see the heavens by naked eye or 
by a telescope, the determination of the falling trend of the 
limiting magnitude, their cross correlation with statistical 
parameters, the determination of the visibility of astronom- 
ical phenomena, the recognition of the areas of a territory 
where the perception of the Universe is more endangered or 
where astronomical observations or popularisation can be 
still acceptably made. In fact (i) the altitude of a site not 
only acts on the levels of sky glow but also has non-negligible 
effects on stellar extinction, (ii) the natural sky brightness 
needs to be accounted for when computing the total sky 
brightness in low polluted sites, (iii) the relation between 
total sky brightness and visual limiting magnitudes is not 
linear, being related to eye capability to see a point source 
towards a bright background. 

Here we extend the method of Paper I in order to be 
able to map naked eye star visibility and telescopic limiting 
magnitudes across large territories. As in Paper I, we eval- 
uate the upward light flux based on DMSP satellite data 
and compute the maps modelling the light pollution prop- 
agation in the atmosphere with Garstang models. They as- 
sume Rayleigh scattering by molecules and Mie scattering by 
aerosols and take into account extinction along light paths 
and Earth curvature. In this paper we take into account the 
altitude of each land area from GTOPO30 elevation data, 
the natural sky brightness in the chosen sky direction, based 
on the Garstang (1989) models, the eye capability or tele- 
scopic limiting magnitudes based on Garstang (2000b) and 
Schaefer (1990) approach, the stellar extinction in the vi- 
sual photometric band based on Snell & Heiser (1968) and 
Garstang (1989) formulae. For near zenith sky directions we 
also take into account mountain screening. 

In section ^| we describe our improvements to the map- 
ping technique. In section^ we deal with input data, describ- 
ing GTOPO30 elevation data, updating the reduction of 
satellite radiance data, summarizing the atmospheric model. 
In section ^ we present the maps of naked eye star visibility 
and total sky brightness in V band in Europe at zenith with 
a resolution of approximately 1 km, we compare map pre- 
dictions with measurements of sky brightness and limiting 
magnitude and we discuss the screening effects. Section H 
contains our conclusions. 



2 MAPPING TECHNIQUE 
2.1 Artificial sky brightness 

The total artificial sky brightness in a given direction of the 
sky in (a/, y') is: 



b(x', y 



e(x, y)f((x, y), {x , y')) dx dy , 



(1) 



where e(x, y) is the upward emission per unit area in (a;, y), 
f((x, y), (a/, y 1 )) is the light pollution propagation function, 
i.e. the artificial sky brightness per unit of upward light emis- 
sion produced by unit area in (x, y) in the site in (x 1 , y'). 
When upward light flux is obtained from satellite measure- 
ments, the territory is divided into land areas with the same 
positions and dimensions as projections on the Earth of the 
pixels of the satellite image and each land area is assumed to 
be a source of light pollution with an upward emission e x , y 



proportional to the radiance measured in the corresponding 
pixel multiplied by the surface area (see eq. 35). The total 
artificial sky brightness at the centre of each area, given by 
the expression (Jlh , becomes: 



(2) 



for each pair and (h, k), which are the positions of the 
observing site and the polluting area on the array. 

In paper I the method of mapping artificial sky bright- 
ness has been applied (i) computing brightness at sea level, 
(ii) assuming sources at sea level and (iii) assuming that the 
upward emission function has the same shape everywhere. 
In this case the light pollution propagation function / de- 
pends only on the distance between the site and the source, 
and on some details which are assumed to be the same ev- 
erywhere, such as the shape of the emission function, the 
atmospheric distribution of molecules and aerosols, their op- 
tical characteristics in the chosen photometric band and the 
direction of the sky observed. Eq. however, it is not 
a convolution because the distance between pairs of points 
depends on the latitude in the used latitude/longitude pro- 
jection. If assumption (i) is relaxed, but assumptions (ii) 
and (iii) are retained, a reasonable computational speed can 
be still obtained evaluating once for each latitude the array 
f(di-h,j-k, hm), where di—hj-k is the distance between the 
site and the polluting areas and m is an index which dis- 
cretizes the altitude h of the site. Both di-h,j-k and h m are 
computed inside a reasonable range (e.g. a circle with 200km 
of radius and the altitude of the highest mountain) . Then, all 
bi t j at nearly the same latitude can be rapidly obtained from 
eq. [^interpolating the array f(di-h,j-k, h m ) at the elevation 
h(i,j) of each site. If assumptions (ii) and (iii) are relaxed it 
becomes necessary to evaluate f((xi,yj),(xh,yk)) for each 
pair of points so that the computations become slower and 
at the moment can be applied only to small territories. For 
this reason, we maintained assumptions (ii) and (iii) when 



computing the maps of Europe in sec. 4.1 



We obtained the propagation functions f(di-h,j-k> h m ) 
or f((xi,yj), (xh,yk)) with models for the light propagation 
in the atmosphere based on the modelling technique intro- 
duced and developed by Garstang (1986, 1987, 1988, 1989a, 
1989b, 1989c, 1991a, 1991b, 1991c, 1992, 1993, 2000a) and 
also applied by Cinzano (2000a, 2000b). The propagation 
function /, expressed as total flux per unit area of the tele- 
scope per unit solid angle per unit total upward light emis- 
sion, is obtained for each set of indexes integrating along the 
line of sight: 

/= /(A»W/m(w)+A(fc)/.(w))(i+l>s)i(V,«)fi(«)d«, (3) 

where /3m(h), /3&(h) are respectively the scattering cross sec- 
tions of molecules and aerosols per unit volume at the el- 
evation h(u) along the line of sight, f m (u), / a (^) are their 
normalized angular scattering functions, £1 (u) is the extinc- 
tion of the light along its path to the telescope and i(ip, s) is 
the direct illuminance per unit flux produced by each source 
on each infinitesimal volume of atmosphere along the line- 
of-sight of the observer. The scattering angle uj, the emission 
angle ip, the distance s of the volume from the source and 
the elevation h of it, depend on the altitudes of the site 
and the source, their distance, the zenith distance and the 
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Figure 1. Screening by terrain elevation. 



azimuth of the line-of-sight, and the distance u along the 
line of sight, through some geometry described for curved 
Earth and nonzero altitude in Garstang (1989, eqs. 4-11). 
The (1 + Ds) in eq. (^) is a correction factor which take into 
account the illuminance due to light already scattered once 
from molecules and aerosols which can be evaluated with the 



approach of Treanor (1973) as extended by Garstang (1984, 
1986, 1989). Details for curved Earth can be found in the 
last paper (eq. 23) as well a discussion about the error in 
neglecting third and higher order scattering which can be 
significant for optical thickness higher than about 0.5. We 
refer the reader to Paper I and the cited papers for details. 

We can account for screening effects setting the illumi- 
nance per unit flux in eq. (0) to: 



(4) 



where there is no screening by Earth curvature or by terrain 
elevation and zero elsewhere. Here I(tp) is the normalized 
emission function giving the relative light flux per unit solid 
angle emitted by each land area at the zenith distance ip, s is 
the distance between the source and the considered infinites- 
imal volume of atmosphere and £2 is the extinction along 
this path. We considered every land area as a point source 
located in its centre except when i = h, j = k in which case 
we used a four points approximation (Abramowitz & Ste- 
gun 1964). Taking into account screening effects requires us 
to check for each point along the line of sight whether the 
source area is screened by terrain elevation or not, taking 
into account Earth curvature. This can be done by deter- 
mining the position of the foot of the vertical of the consid- 
ered point and then computing for every land area crossed 
by a line connecting this foot and the source area, the quan- 
tity cot tfj, where if) is defined as in figure [l] and depends on 
the elevation A of the land area, the distance D of its center 
from the center of the source area and the Earth radius E: 



cot ift = 



(A + E) - (h + E) cos(D/E) 



(h + E) sin{D/E) 



(5) 



Then we can determine X 
screening elevation h s : 



max(— cot ip) and from it the 



A + E 

cos(7> IE) - X sin(D* / E) 



(6) 



where D* is the distance between the source area and the 
foot of the vertical, and h s is computed over the sea level. 
The illuminance i in eq. is set to zero when the elevation 
of the considered point is lower than the screening elevation. 
For lines of sight pointing toward the zenith the evaluation 
of the screening elevation can be done once for each pair site- 
source. However, even in this faster case the computational 
time required by the screening evaluation for each source 
area around each site is huge, so we accounted for screening 
effects due to terrain elevation only in the small maps of 
section 4.4 accounting in the other maps only for screening 
by Earth curvature as described by Garstang (1989, eq. 12- 
13). 



2.2 Natural sky brightness 

The mapping of the total sky brightness requires the eval- 
uation of the natural sky brightness under the atmosphere, 
i.e. as actually observed from the ground, in the direction of 
the line-of-sight. 

We assumed as Roach & Meinel (1955) and Garstang 
(1989) that natural sky brightness is produced by (i) light 
from a layer at infinity due mainly to integrated star light, 
diffused galactic light and zodiacal light, and (ii) light due to 
airglow emission from a van Rhijn layer at height of 130 km 
above the ground. The first, b s , depends on the equatorial co- 
ordinates of the observation point and, for the zodiacal light, 
on the time t. The second depends on the angle at which the 
layer is observed and on the layer brightness, b VI , which in 
turn depends on some factors like the geographical latitude 
L, the solar activity S in the previous day, and the time 
T after twilight. Extending the Garstang (1989) approach, 
we assumed that the natural sky brightness bout outside the 
scattering and absorbing layers of the atmosphere, at zenith 
distance z and azimuth u, is given by: 



&out(z, u,L,T, S) = b B (a,8,t) + 



b VI (L,T,S) 
(1 -0.96 sin 2 z)V2 



(7) 



where a is the right ascension and 5 the declination of the 
observed area of the sky which depends on the zenith dis- 
tance 2 and azimuth u) through the observation time t. From 
Walker (1988) we know that the sky brightness decays to a 
nearly constant level after some hours from the astronomical 
twilight due probably to the recombination of ions excited 
during the day by the solar radiation, so we will refer the 
night brightness always to some hours after the twilight in 
order that the dependence on T disappear. The dependence 
of & vr on the solar activity can be expressed in a very rough 
approximation as Cinzano (2000c), based on the measure- 
ments of Walker (1988), Cannon 1987 in Krisciunas et al. 
(1987), Krisciunas (1990, 1997): 



: &vr,o(L)-10' 



C cos(2iT t —^L) 



(8) 



where & vr ,o is the value at mean solar activity, P is the aver- 
age period of the solar cycle, to is the epoch of a maximum, 
t is the time from the epoch, and C is a constant. A good 
correlation was found by Walker (1988, see also Krisciunas 
1997, eq. 5) between the sky brightness and the observed 
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10.7 cm solar radio flux for the day preceding the obser- 
vations. The value of b s in a given direction on the sky, 
depending on the observation time, can be significant when 
maps are made to quantify the visibility of an astronomical 
phenomenon. When the purpose is to evaluate the generic 
capability to see the stars, we can assume b s constant and 
given by its average value at the considered latitude. In this 
case 6 out does not depend on the azimuth. 

The natural sky brightness 6 na t observed at altitude A 
over sea level is given by the sum of (i) the light directly com- 
ing to the observer fed, (h) the light scattered by molecules 
& m and (iii) the light scattered by aerosols fe a : 



fenat = fed + fern + ba 



(9) 



We computed these components using the model presented 
by Garstang (1989). 

The directly transmitted light arrives at the observer 
after extinction along the line of sight (Garstang 1989): 



fed = feout£3 



(10) 



where £3 is the atmospheric extinction in the path from in- 
finity to the observer. 

The light scattered by molecu les by Rayleigh scattering 
for the atmospheric model of sec. 3.3 is (Garstang 1989): 



3(1 + GQfeoutfli 
16 exp(cH) 



i s (u, 2) exp( — ch)^4,du (11) 



is{u, z) — / (2 + 2/i 2 cos 2 2 + (1 — fi 2 ) sin 2 z) f(fi)dfi (12) 

^ + ^L(0.04 + 0.96 M 2 )- 1/2 )e 3J D S 2 (13) 

'out Oout 

1 + ll.lLRT/3 m exp(-aff)a~ V" 1 exp(-ah) (14) 



/O) 

D s2 



where u is the integration variable along the line of sight, 
/3 m is the scattering cross section of molecules per unit vol- 
ume, G is the ground reflectivity (assumed 15%), c is the 
molecular inverse scale height, a is the aerosol inverse scale 
height, H is the elevation of the land, h is the elevation of 
the infinitesimal volume at u, A is the elevation of the site, 
H is an integration variable, K is the Garstang atmospheric 
clarity coefficient, £4 is the atmospheric extinction in the 
path from the scattering volume to the observer. 

The ligh t sc attered by aerosols for the atmospheric 
model of sec. 



3.3 is (Garstang 1989): 



fea = feoutd aJ 4M£3.Ds3 (15) 

d a = ll.ll-ft:/3 m exp(--cii')a~ 1 exp(-a J 4) (16) 
D s3 = l + 0.5d a A M (17) 

where Am is the airmass given by eq. (|iil). 

Total sky brightness is fe-rij = bij + fo na t. We expressed 
it in photon radiance in ph cm~ 2 s _1 sr _1 or in magnitudes 
per arcsec 2 with the Garstang (1986, 1989) relation: 



Vi. 



41.438 - 2.51og& Ti j 



(18) 



We evaluated 6 S and 6 vr by fitting the predictions for the nat- 
ural sky brightness at zenith with specific observations made 
in unpolluted sites after some hours from sunset and reduced 
fevr approximately with eq. (^) to average solar activity. We 
assumed 6 a and 6 vr ,o constant everywhere in Europe. 



2.3 Naked eye star visibility 

We obtained the map of naked eye limiting magnitude as 
an array rriij, giving the limiting magnitude at each grid 
point, from the array Vij, giving the total sky brightness, 
as described. 

The illumination i in lux perceived by the eye from 
a source which is at the threshold of visibility to an ob- 
server when the brightness of observed background is bobs 
in nanolambert and the stimulus size, i.e. the seeing disk di- 
ameter, is 8 in arc minutes, was given by Garstang (2000b) 
based on measurements of Blackwell (1946) and Knoll, Tou- 
sey and Hulburt(1946): 

h = c 1 (l + fc 1 fe 1/2 ) 2 (l + a 1 6» 2 +y 1 6^ s '' 2 ' 



obs 

ia = c 2 (l + k 2 b 1/2 ) 2 {l + a 2 6 2 +y 2 b z o l s 6 
i = iii 2 /(ii+i 2 ) 



2\ 



with a = 3.451 x 10 



-9 



1.51 x 10" 



yi 



c 2 = 4.276 x 10" 
2.0 x 10" 5 , y 2 = 



(19) 
(20) 
(21) 
0.109, 



, fei 
1.29 x 10" 



2i = 0.174, 22 = 0.0587, ai = 2.35 x 10"" 4 , a 2 = 5.81 x 10 -3 . 
The last equation is an artifact introduced by Garstang in 
order to put together smoothly the two components ii and 
i 2 , related respectively to the thresholds of scotopic and pho- 
topic vision, obtaining the best fit with cited measurements. 

The observed background b b s in eqs. (19-20) is related 
to the night sky background b V i S in the visual band from 
(Garstang 2000b): 



b vis /(F,F sc F c] 



(22) 



where F a takes into account ratio between average pupil 
area of the Knoll, Tousey, Hulburt and Blackwell observers 
and the pupil area of the assumed observer, Fsc takes into 
account the Stiles- Crawford effect, F c ^ allows for the differ- 
ence in colour between the laboratory sources used in deter- 
mining the relationships between i and b and the night sky 
background. Sky brightness b v j s in visual band expressed in 
nanolambert can be obtained from sky brightness Vij in V 
band, expressed in mag/ arcsec 2 , inverting Garstang (1986, 
1989) relation: 

b vis = io-°- 4 ^- 26 - 346 > (23) 

The illumination i 1 produced over the atmosphere by a 
star at the threshold visibility are related to the threshold 
illuminations ii,i 2 obtained from eqs. (19-20) for scotopic 
and photopic vision from (Garstang 2000b): 

i' = iW 2 /(i'i+i 2 ) (24) 
ii = F & ,iFsc,iFcB,iFe,iF s ,iii (25) 
i' 2 = F,, 2 F S c,2F CB , 2 F c , 2 F B!2 i 2 (26) 

where F a and Fsc are defined as before, F cs allows for the 
difference in color between the laboratory sources and the 
observed star, F c allows for star light extinction in the terres- 
trial atmosphere, taking into account that star magnitudes 
are given outside the atmosphere, F s allows for the acuity 
of any particular observer, defined so that F s < 1 leads 
to a lower threshold i and therefore implies an eye sensi- 
tivity higher than average due possibly to above average 
retinal sensitivity, observing experience or an above average 
eye pupil size. The illumination i' expressed in lux can be 
converted into magnitudes (Allen 1973 p. 197): 



-13.98 - 2.5 log % 



(27) 
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If po is the pupil diameter used by the average of the 
Knoll, Tousey, Hulburt and Blackwell observers, who are as- 
sumed to have been age Ao — 23, and p is the pupil diameter 
of an observer aged A when the sky background is b b s , then 
F a = pl/p 2 from eq. (6) of Garstang (2000b) becomes: 



F 



0.534-0.00211Aq -(0.236-0.00127Ao)g 
0.534-0.0021L4-(0.236-0.00127A)q 
q = tanh(0.40 log fo obs - 2.20) 



(28) 
(29) 



The Stiles-Crawford effect, due to the decreasing efficiency 
in detecting photons with the distance from the center of 
the pupil, produces a non linear increase of sensibility when 
the eye pupil increases and can be taken into account with 
equations from Schaefer (1990) modified as pointed out by 
Garstang (2000b) and Schaefer (1993). We neglected this 
effect which must be accounted for telescopic observations 
or whenever a larger precision is needed. 

Differences in colour between the eye sensitivity curve 
and photometer sensitivity curve used in determining i can 
be corrected with (Schaefer 1990): 



F cs ,i 
F cs , 2 



10 _0.4(l-(S-V)/2) 
1 



(30) 
(31) 



where we assumed as a typical star color index B— V = 
0.7. Differences in colour between the laboratory background 
and the night sky background F c b can be corrected with the 
same formula. The colour of night sky is difficult to evalu- 
ate when there is light pollution. Cinzano & Stagni (2000) 
showed that the sky becomes redder far from sources. How- 
ever near predominant sources, like large cities, where the 
extinction is negligible and aerosol scattering is large, the 
colour index is related to the colour of the integrated lamp 
spectra. We assumed here B— V = 0.7 on average but when 
emission spectra of each land area become available, it will 
be possible to obtain the colour index of the night sky point 
by point computing maps of total sky brightness separately 
for B and V bands like in Paper I. Stellar extinction in the 
atmosphere F c is computed from eq. ( p^ ) and from the V 
band vertical extinction (sec. 3.3) corrected approximately 



for the night vision as Schaefer (1990). In the computation 
of i[ and i' 2 must be used respectively the correction fac- 
tors for scotopic and photopic vision. The reader is referred 
to Schaefer (1990) and Garstang (2000b) for an extensive 
discussion. 



2.4 Telescopic limiting magnitude 

We can also obtain maps of telescopic limiting magnitude 
for a given instrumental setup. This could be useful for am- 
ateurs observational campaigns. In this case we must replace 
the image size 6 by M9 in eq. (|l9|), where M is the magni- 
fication of the telescope. 

The observed background 6 b s is related to the night 
sky background under the atmosphere 6 v is from (Garstang 
2000b): 



= 6 v is/ (F h F t F p F a .F nl FscF c 



(32) 



where F & takes into account the ratio of the area of the tele- 
scope to that of the naked eye, Fsc takes into account the 
Stiles-Crawford effect, _F c b allows for the difference in colour 



between the laboratory sources used in determining the re- 
lationships between i and b and the night sky background, 
Fb takes into account that one eye is used in telescopic ob- 
servations, while binocular vision was used in obtaining the 
relations between i and b, Ft allows for the loss of light in 
the telescope, Ft being the reciprocal of the transmission t 
through the telescope and eyepieces, F p allows for the loss of 
light if the telescope exit pupil is larger than the eye pupil, 
F m allows for the reduction of the sky brightness by the 
telescope magnification 

The illuminations perceived from a star are related to 
the illuminations given by eqs. (19-20) for scotopic and pho- 
topic vision from (Garstang 2000b): 



ii = F a ,iFsc,iF CSl iFe,iF s ,iF M Ft,iFp,rii 

i'2 = F a ,2Fsc,2F CSi 2F e ,2F Si 2F b ,2Ft,2F Pi 2«2 



(33) 
(34) 



where F cs , F e and F s has been already discussed in sec. 2.3 



We refer the reader to Schaefer (1990) and Garstang 
(2000b) for the formulae and further discussions. Note that, 
as pointed out by the last author, Schaefer's F r is not needed 
because the image size was already included in eq. (p"s|). 



3 INPUT DATA 
3.1 Altitude data 

As input elevation data we used GTOPO30, a global digi- 
tal elevation model (DEM) by the U.S. Geological Survey's 
EROS Data Center. Details have been given by Gesch et al. 
(1999). This global data set covers the full extent of latitude 
and longitude with an horizontal grid spacing of 30-arc sec- 
onds as does our composite satellite image. From the global 
16-bit DEM (21,600 rows by 43,200 columns), provided as 
16-bit signed integer data in a simple binary raster, we cut 
an array of 4800x4800 pixels covering the same area as our 
satellite image. The vertical units represent elevation in me- 
ters above mean sea level which ranges from -407 to 8,752 
meters. We reassigned a value of zero to ocean areas, masked 
as "no data" with a value of -9999, and to altitudes under 
sea level. 

GTOPO30 is based on data derived from 8 sources of 
elevation information, including vector and raster data sets: 
(i) Digital Terrain Elevation Data (DTED) a raster topo- 
graphic data base with a horizontal grid spacing of 3-arc 
seconds (approximately 90 meters) produced by the Na- 
tional Imagery and Mapping Agency (NIMA); (ii) Digital 
Chart of the World, a vector cartographic data set based on 
the l:l,000,000-scale Operational Navigation Chart series 
products of NIMA; (iii) USGS 1-degree DEM's with an hor- 
izontal grid spacing of 3-arc seconds; (iv) Army Map Service 
l:l,000,000-scale paper maps (AMS) digitized by Geograph- 
ical Survey Institute (GSI) of Japan; (v) International Map 
of the World l:l,000,000-scale (IMW) digitized by GSI for 
the Amazon basin; (vi) digitized Peru l:l,000,000-scale map 
to fill gaps in source data for South America; (vii) Manaaki 
Whenua Landcare Research DEM with a 500-meter hori- 
zontal grid spacing for New Zealand; (viii) Antarctic Digital 
Database (ADD) under the auspices of the Scientific Com- 
mittee on Antarctic Research. 

The absolute vertical accuracy varies by location ac- 
cording to the source data and at the 90 percent confidence 
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level is 30 m for DTED, 160 m for DCW, 30 m for USGS 
DEM, 250 m for AMS maps, 50 m for IMW maps, 500 m for 
Peru map, 15 m for N.Z. DEM and not available for ADD 
(Gesch et al. 1999). For many areas the relative accuracy is 
probably better than the estimated absolute accuracy. 

As discussed by Gesch et al. (1999), due to the nature of 
the raster structure of the DEM, small islands in the ocean 
less than approximately 1 square kilometre are not repre- 
sented. Nonetheless the error in assuming them at sea level 
is usually small because their limited size do not allow very 
high elevations. Not all topographic features that one would 
expect to be resolved at 30-arc second grid spacing are rep- 
resented but this grid spacing is appropriate for the areas 
derived from higher resolution DEM's. Changes in detail of 
topographic information are evident at the boundary be- 
tween two sources, even if the mosaicing techniques smooth 
the transition areas. Artefacts due to the production method 
are plainly visible in some areas even if their magnitudes in 
a local area are usually well within the estimated accuracy 
for the source. Some production artefacts are also present 
in the areas derived from the vector sources. Small artificial 
mounds and depressions may be present in localized areas, 
particularly where steep topography is adjacent to relatively 
level areas, and the data were sparse. Additionally, a "stair 
step" (or terracing) effect may be seen in profiles of some 
areas, where the transition between contour line elevations 
does not slope constantly across the area but instead is cov- 
ered by a flat area with sharper changes in slope at the 
locations of the contour lines. 



3.2 Upward flux data 

Upward flux data have been obtained from the Operational 
Linescan System (OLS) carried by the U.S. Air Force De- 
fense Meteorological Satellite Program (DMSP) satellites. 
This is an oscillating scan radiometer with low-light visible 
and thermal infrared (TIR) imaging capabilities. At night 
the OLS, carrying a 20 cm reflector telescope, uses a Photo 
Multiplier Tube (PMT) to intensify the visible band sig- 
nal which have a broad spectral response covering the range 
for primary emissions from the most widely used lamps for 
extern al lighting. Details are described in Paper I, Lieske 
( [L98l| ), Elvidge et al. (1997a, 1997b, 1997c, 1999). 

The collection of special DMSP data, used in Paper 
I and here to assemble a cloud-free composite image cal- 
ibrated to top-of-atmosphere radiances, has been obtained 
after a special requests to the Air Force made by the U.S. De- 
partment of Commerce, NOAA National Geophysical Data 
Center (NGDC), which serves as the archive for the DMSP 
and develops night time lights processing algorithms and 
products. The primary reduction steps include (see Paper I 
and Elvidge et al. 1999): 

1) acquisition of special OLS-PMT data at a number of re- 
duced gain settings (24dB, 40dB, 50dB) to avoid saturation 
on major urban centres and, in the same time, overcome 
PMT dynamic range limitation (Our request was granted 
for the darkest nights of lunar cycles in March 1996 and 
January- February 1997.). On board algorithms which ad- 
just the visible band gain were disabled. 

2) establishment of a reference grid with finer spatial reso- 
lution than the input imagery; 

3) identification of the cloud free section of each orbit based 



on OLS-TIR data; 

4) identification of lights, removal of noise and solar glare, 
cleaning of defective scan lines and cosmic rays; 

5) projection of the lights from cloud-free areas from each 
orbit into the reference grid; 

6) calibration to radiance units using prior to launch calibra- 
tion of digital numbers for given input telescope illuminance 
and gain settings in simulated space conditions; 

7) tallying of the total number of light detections in each 
grid cell and calculation of the average radiance value; 

8) filtering images based on frequency of detection to remove 
ephemeral events; 

9) transformation into latitude/longitude projection with 
30"x30" pixel size; 

10) cutting of the requested portion of the final image 
(we used an image of 4800x4800 pixel corresponding to 
40° x 40° starting approximately at longitude 10°30' west and 
latitude 72°north). 

We improved map predictions by applying to the com- 
posite satellite image a mild deconvolution with the Lucy- 
Richardson algorithm. In fact, (i) the effective instantaneous 
field of view is larger (2.2 km at nadir to 5.4 km at the scan 
edges) than pixel-to-pixel ground sample distance (GSD) 
maintained by the along-track OLS sinusoidal scan and the 
electronic sampling of the signal from the individual scan 
lines (0.56 km), (ii) most of the data received by NGDC 
has been "smoothed" by on-board averaging of 5 x 5 pixel 
blocks, yielding data with a GSD of 2.8 km, and (iii) data 
of more orbits have been tallied. A mild deconvolution al- 
lows partial recovery of the smeared radiance and allows bet- 
ter predictions for sites near strong sources like cities where 
spreading in distribution of upward emission could have an 
effect on map results. The point spread function has been 
obtained searching for isolated nearly-point sources and the 
deconvolution has been applied to smaller subimages. We 
plan in future analysis to download from the satellite the 
original high-resolution data in order to properly deconvolve 
single orbit data before tallying. We also plan in future re- 
ductions to correct data for atmospheric extinction before 
tallying and to use only data from areas not very far from 
nadir in order to avoid effects of the shape of the upward 
emission function when changing the observation angle. As 
showed in Paper I, however, this is only a second order ef- 
fect due to opposite contributions of extinction and emission 
function shape. 

Calibrated upward flux measurements can be obtained 
based on pre-fly irradiance calibration of OLS PMT as 
described in Paper I. If r is the energetic radiance in 
[l0 10 Wcm" 2 sr _1 ] measured by the OLS-PMT, the upward 
light flux e in [V band photons s _1 ] is given by eq. (28), 
(29) and (30) of Paper I: 



f °°T v h\d\ < A > 10°- 4Am cos(0 
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where 7(?/>) is the normalized upward emission function as 
defined in paper I, I the latitude of the land area, 7?t the 
average Earth radius in km and Ax the pixel size in arcsec, 
Tv and Tpmt the sensitivity curves respectively of V band 
and PMT detector, 7a the energy spectrum of the emission 
from the chosen land area, h the Planck constant, c the ve- 
locity of light and < A > the effective wavelength of the 
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combination of the sensitivity curves of the PMT and the 
calibration source. We assumed an average vertical extinc- 
tion Am = 0.33 mag V for all measured land areas ne- 
glecting their elevation and solved the integrals in eq. ( |35[ ) 
constructing, as in Paper I, an approximate synthetic spec- 
trum for a typical night-time lighting roughly assuming that 
50 per cent of the total emitted power be produced by High 
Pressure Sodium (HPS) lamps (SON standard) and 50 per 
cent by Hg vapour lamps (HQL). 

We assumed that all land areas have on average the 
same normalized emission f unctio n, given by the parametric 
representation of Garstang (1986) in eq. (15) of Paper I. This 



has been tested by Garstang and by Cinzano (2000a) with 
many comparisons between model predictions and measure- 
ments and in Paper I by studying in a single orbit satel- 
lite image the relation between the upward flux per unit 
solid angle per inhabitant of a large number of cities and 
their distance from the satellite nadir, which is related to 
the emission angle ip. See Paper I for a detailed discussion. 



3.3 Atmospheric data and stellar extinction 

In order to compute extinctions along the light paths and 
scatterings, we need atmospheric data or an atmospheric 
model for the considered territory. In principle what we 
need is a set of functions giving, for each triplet of lon- 
gitude, latitude and elevation (x,y,h), the molecular and 
aerosol scattering coefficients per unit volume of atmosphere 
Pm(x, y, h) and I3 a (x,y,h), and the aerosol angular scatter- 
ing function f a (uj,x,y,h). The molecular angular scattering 
function fm(u>) is known because it is Rayleigh scattering. If 
discretised in arrays, they should preferably have the same 
grid spacing of our upward flux and elevation data. The 
atmospheric data or model would need to refer to condi- 
tions for typical clean nights at every point and contain any 
other information on denser aerosol layers, volcanic dust and 
the Ozone layer. At the moment this is not at our disposal. 
Moreover, applying typical condition in every land area we 
risk mixing effects due to light pollution with effects due to 
gradients of atmospheric conditions in typical nights. So we 
applied the same standard atmospheric model everywhere, 
neglecting geographical gradients and local particularities as 
in Paper I. Here we resume the simple atmospheric model. 

1) We assumed the molecular atmosphere in hydr ostati c 
equilibrium under the gravitational force as Garstang ( 1986 ) 



with an inverse scale height c = 0.104 km - , a molecular 
density at sea level iV m ,o = 2.55 x 10 19 cm -3 and a constant 
integrated Rayleigh scattering cross section in V band a m = 
4.6 x 10 -27 cm 2 molecule -1 . The scattering cross section 
per unit volume is fi m ,o = ^m.o^n. Note that in sec. 4.2 of 
Paper I the letters "B" and "V" of the photometric bands 
of the Rayleigh cross sections have been exchanged and that 
the units are [cm 2 molecule -1 ] . 

2) We assumed the atmospheric haze aerosols num- 
ber densit y dec reasing exponentially with the altitude as 
Garstang (1986) and Joseph et al. (1991) neglecting the 



presence of sporadic denser aerosol layers, volcanic dust and 
the Ozone laye r, studied by Garstang (1991a, 1991c). As 
Garstang (1986), the inverse scale height of aerosols was as- 
sumed to be a — 0.657 + 0.059-Zf. Aerosol content was given 
using the Garstang atmospheric clarity parameter K which 



measures the relative importance of aerosol and molecules 
for scattering light in V band: 



K 



An.oll.lle-oH 



(36) 



where H is the altitude of the ground level above sea level. 
The typical normalized angular scattering function for at- 
mospheric haze aerosols was assumed to be g iven by the 
function tabulated by Mc Clatchey et al. (197£) as interpo- 
lated by Garstang (1991a). 



The stellar extinction at zenith in magnitudes for a 
site at altitude A is given for this atmospheric model by 
Garstang (1986): 



Am = 1.0857/3 m , e" 
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(37) 



The stellar extinction at zenith distance z can be obtained 
from the air mass ^4m given by Snell & Heiser (1968): 

Am = sec z ■ 



g(sec z 



(38) 



with g — 0.010 as chosen by Garstang (1986) to reproduce 
the table I, column 3 of Allen (1973) to better than 0.1 air 
masses at zenith distance 85°. The atmospheric extinctions 
£i, £2, £3, £4 of sec. |^ are given for this atmospheric model 
by Garstang (1989, eqs. 18-22). We can associate the at- 
mospheric conditions with other observable quantities like 
the horizontal visibility (Garstang 1989, eq. 38), the optical 
thickness r (Garstang 1986, eq. 22) and the Linke turbidity 
factor for total solar radiation (Garstang 1988). 



4 RESULTS 

4.1 Maps of total sky brightness and star visibility 

We present as an application the map of zenith night sky 
brightness and naked eye star visibility in Europe. 

Fig. ^ shows the total sky brightness at the zenith in 
V photometric astronomical band (Johnson 1955). Colour 
levels from brown to white correspond to total sky bright- 
ness of: <17.5, 17.5-18, 18-18.5, 18.5-19, 19-19.5, 19.5-20, 
20-20.5, 20.5-21, 21-21.5, >21.5 V mag/arcsec 2 . The map 
was computed for clean atmosphere with aerosol clarity 
K = 1, corresponding to a vertical extinction in V band 
of Am = 0.33 mag at sea level, Am = 0.21 mag at 1000m 
o.s.l., Am = 0.15 mag at 2000m o.s.l., horizontal visibility 
at sea level Ax — 26 km, optical depth r = 0.3 so double 
scattering approximation is adequate. Each pixel is 30" x30" 
in size in longitude/latitude projection. Computation of the 
natural sky brightness is based on measurements at Isola 
del Giglio (Italy), a quite dark site according to the World 
Atlas of Artificial Sky Brightness (Cinzano et al. in prep.), 
giving V — 21.74 ± 0.06 mag arcsec -2 in V band in 1999 
at 200m o.s.l.. At that time average solar activity was re- 
ported. The map was rescaled from 1996-1997 to 1998-1999 
adding Am = 0.28 mag/arcsec 2 to the total sky brightness 
obtained from OLS-PMT pre-fly radiance calibration. This 
correction was obtained fitting a straight line y = Am + x 
to the observed versus predicted data points of sec. 4.2. The 
difference Am = 0.28 ± 0.10 mag/arcsec 2 is possibly due 
to the growth of light pollution and agree within errorbars 
with Cinzano (2000c). 
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Figure 2. Total night sky brightness in Europe in V band for aerosol content parameter K = 1. 
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Figure 3. Naked eye limiting magnitude in Europe for aerosol content parameter K = 1. 
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The limiting magnitude is a statistical concept (Black- 
well 1946; Schaefer 1990). A number of random factors affect 
eye measurements like the individual eye sensitivity, the dif- 
ference from the average eye pupil size, the capability to use 
averted vision, the experience which make an observer con- 
fident of a detection at a probability level different from the 
others (Schaefer 1990 reports a difference of 1 magnitude 
from 10% to 90% detection probability, corresponding re- 
spectively to the fainter suspected star and the fainter surely 
visible star), the lenghts of time for which the field has been 
observed (Schaefer 1990 reports roughly half a magnitude 
from 6-seconds to 60-seconds observation) . Fig. ^ shows the 
centre of the gaussian distribution of the naked eye limiting 
magnitude in Europe at the zenith obtained from the map 
of fig. |^. Colour levels from pink to black correspond to: 
<3.75, 3.75-4.00, 4.00-4.25, 4.25-4.50, 4.50-4.75, 4.75-5.00, 
5.00-5.25, 5.25-5.50, 5.50-5.75, 5.75-6.00, >6.00 V magni- 
tudes. Limiting magnitudes are computed for observers of 
average experience and capability F s = 1, aged 40 years, 
with the eyes adapted to the dark, observing with both 
eyes. Observer experience can be accounted with eq. (20) 
of Schaefer (1990). 

Original maps are 4800x4800 pixel images saved in 16- 
bit standard fits format with fitsio Fortran- 77 routines de- 
veloped by HEASARC at the NASA/GSFC. They have been 
analysed with FTOOLS 4.2 analysis package by HEASARC 
and with Quantum Image 3.6 by Aragon Systems. Read- 
ers should consider the distinction between grid spacing and 
resolution. The resolution of the maps, depending on re- 
sults from an integration over a large zone, is greater than 
resolution of the original deconvolved images and is gener- 
ally of the order of the distance between two pixel centres 
(30"x30", i.e. less than 1 km). Country boundaries are ap- 
proximate. 

4.2 Comparison with total sky brightness 
measurements 

Fig. ^ shows a comparison between predictions of total night 
sky brightness and measurements in Europe in V band in 
1998 and in 1999. A detailed comparison requires measure- 
ments taken (i) at a large number of sites, (ii) in nights 
with the same vertical extinction and horizontal visibility as- 
sumed in the map computation, (iii) in many similar nights 
in order to smooth atmospheric fluctuations by averaging, 
(iv) under the atmosphere, i.e. as actually observed from 
the ground without any extinction correction applied, (v) 
in the same period in which the satellite image was taken 
in order to minimize uncertainties given by the fast growth 
rate of artificial sky brightness, (vi) in the same photomet- 
ric band for which maps are computed, (vii) with accurate 
geographical positions. As in Paper I, due to the scarcity 
of measurements of sky brightness associated to measure- 
ments of extinction, or to any other index of the atmospheric 
aerosol content, we used all available measurements taken in 
clean or photometric nights even if extinction was not avail- 
able or not exactly the required one (Catanzaro & Cata- 
lano 2000; Cinzano 2000a; Delia Prugna 1999; Favero et al. 
2000; Piersimoni, Di Paolantonio & Brocato 2000; Poretti 
& Scardia 2000; Zitelli 2000, Falchi 1998). Differently from 
Paper I we did not need to subtract an assumed natural sky 
brightness from measurements to obtain artificial sky bright- 




predicted brightness 



Figure 4. Measurements of total sky brightness versus map pre- 
dictions in V band. The straight line is the linear regression. 

ness. Errorbars indicate measurement errors which are much 
smaller than the uncertainties produced by fluctuations in 
atmospheric conditions which are unknown. Shifts in mea- 
surements obtained with different instrumental setups could 
also arise, as pointed out in Paper I. The best fit of a straight 
line y — a + bx to the data points (solid line in fig. ^), as- 
suming unknown uncertainties, gives b = 0.99 ± 0.08 and 
a = —0.21 ± 1.58 mag/arcsec 2 . The sigma derived from the 
chi-square assuming that our predictions fit well, a = 0.35 
mag/arcsec 2 , give an estimate of the uncertainty of our pre- 
dictions at a site. When a large number of measurements 
of sky brightness together with their stellar extinction will 
be available, a more precise evaluation of the uncertainty 
become possible. A worldwide CCD measurement campaign 
of both sky brightness and stellar extinction has been orga- 
nized by the International Dark-Sky Association (Cinzano 
& Falchi 2000). 



4.3 Comparison with naked eye limiting 
magnitude measurements 

We compared our predictions of naked eye limiting magni- 
tude with observations taken in Europe in the period 1996- 
1999. Data have been obtained from preliminary results of 
the measurement campaigns set up by a number of organi- 
zations: (i) Dataset A: Operation 'Atlas 1996', Comite Na- 
tional Pour la Protection du Ciel Nocturne, France (Corp 
1998); (ii) Dataset B: 'Gli studenti fanno vedere le stelle', 
Ministero della Pubblica Istruzione - Unione Astrofili Ital- 
iani - Legambiente, Italy (Corbo 2000a, b); (iii) Dataset 
C: Astronomy On-line 'Light Pollution Project', European 
Southern Observatory (Haenel 1999); (iv) Dataset D: 'CCD 
Amateurs Measurements of Night Sky Brightness', Interna- 
tional Dark-Sky Association (IDA) - Italian Section (Falchi 
priv. comm.). A detailed comparison between map predic- 
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Figure 5. Distribution of the O-C residuals for naked eye limiting 
magnitudes in dataset A (upper panel), datasets B and C (middle 
panel), dataset D (lower panel). Best fitting gaussians are also 
shown in the two upper panels. 



tions of naked eye limiting magnitudes and visual estimates 
requires observations made (i) at a large number of sites, 
(ii) by a large number of observers in each site in order to 
have a statistical treatment of eye capabilities, (iii) in nights 
with the same vertical extinction and horizontal visibility as- 
sumed in the map computation, (iv) in many similar nights 
in order to smooth atmospheric fluctuations by averaging, 
(v) in the same period in which the satellite image was taken, 
(vii) with accurate geographical positions (better than 15"). 
The number of measurements from each site at our disposal 
was too little to allow a statistical analysis site by site. More- 
over, excluding the IDA data, the measurements have been 
taken without contemporary photometrical measurements of 
extinction so that the effects of random atmospheric content 
add further uncertainty. Uncertainties in geographical posi- 
tion could also exist. Systematic errors could also arise if the 
majority of observers detected the faint suspected star and 
not the fainter surely seen. We excluded measurements for 
which the observer reported the presence of the moon, un- 




4 5 6 

predicted limiting mag. 

Figure 6. Measurements of naked eye limiting magnitude versus 
map predictions for dataset A (dots), datasets B and C (trian- 
gles) and dataset D (open diamonds). The dotted line separates 
approximately the areas with prominent photopic vision (left) and 
scotopic vision (right). Solid line shows O-C=0 and the dashed 
lines show deviations of ±1 mag. 

clean sky, large light installations in the nearby or for which 
we were unable to determinate the geographical position. 

Fig. H shows the distribution of the observed minus cal- 
culated limiting magnitudes. The upper panel shows dataset 
A (180 sites) made for 50% by active amateurs, for 25% 
by very experienced amateurs and for 25% by beginners. A 
gaussian give a very good fit with R 2 — 0.986, a shift of 
the centre of only Ax = —0.22 ± 0.03 mag and a dispersion 
a = 0.79 ± 0.04 mag. Schaefer (1990) compared 314 visual 
observations with his model of limiting magnitude obtaining 
a nearly-gaussian distribution of errors with an HWHM of 
0.75 mag and a shift of the centre of -0.24 mag. The middle 
panel show measurements coming from datasets B and C (22 
sites), two campaigns devoted to schools and unexperienced 
observers. A gaussian fits well with R 2 = 0.979. The shift of 
its centre toward more luminous stars, Aa; = —0.53 ± 0.03 
mag, and the smaller dispersion a — 0.50 ± 0.03 mag could 
be due to the more homogeneous kind of observers, mainly 
equally unexperienced. The lower panel shows few observa- 
tions from dataset D, a project devoted to advanced ama- 
teurs. Measurements have been obtained by two experienced 
observers searching for the faintest suspected star, with the 
same atmospheric conditions, stated by contemporary mea- 
surements of stellar extinction, and precise geographic po- 
sitions. They show a small standard deviation (0.20 ± 0.07 
mag) and a mean shifted of Aa; « 0.63 mag toward fainter 
visible stars in reasonable agreement with Schaefer (1990) 
formula for expert observers and a 10% threshold. 

Fig. [] shows the observed versus calculated limiting 
magnitudes for dataset A (dots), datasets B and C (tri- 
angles) and dataset D (open diamonds). The dotted line 
separates approximately the areas with prominent photopic 
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vision (left) and prominent scotopic vision (right). Solid line 
correspond to O-C=0 and the dashed lines show deviations 
of ±1 mag. Experienced observers of dataset D seems to 
gain a magnitude going from photopic to scotopic visibil- 
ity. The scotopic observations from the other datasets show 
a larger scatter toward higher magnitudes than photopic 
observations. If confirmed this could be due to the greater 
difficulty to observe in scotopic conditions (e.g. the eye sen- 
sitivity changes with distance from the center of the retina 
and it is very low at the center, see Clark 1990) . Further ob- 
servations made in a selected number of sites with different 
sky brightnesses by a large number of observers in each of 
them could allow a deeper statistical analysis. 



4.4 Screening 

In paper I we neglected the presence of mountains which 
might shield the light emitted from the sources from a frac- 
tion of the atmospheric particles along the line-of-sight of 
the observer. Given the vertical extent of the atmosphere in 
respect to the height of the mountains, the shielding is not 
negligible only when the source is very near the mountain 
and both are quite far from the site (Garstang 1989a; see 
also Cinzano 2000b). However when taking into account al- 
titudes in map computation other two effects of screening 
by terrain elevation can result in plainly visible artefacts, 
even if mainly within the accuracy estimated for the maps. 
In valleys surrounded by mountains, their lower elevations 
relative to nearby mountain edges result, when screening is 
neglected, in a greater sky brightness due to the illumination 
of the part of the line of sight going from valley altitude to 
mountain edge altitude. When proper screening from sur- 
rounding mountains is taken into account, this part of the 
line of sight is shielded and does not contribute to sky bright- 
ness. Due to the extinction along the longer light path, the 
sky brightness is lower than at the mountains edges. An- 
other effect is produced by small elevations of terrain which 
can enhance efficiency of the Earth curvature in shielding 
far sources. This can be evident, e.g. in case of dark areas 
located at a hundred kilometres from very lighted areas. 

Fig. (?] shows the effects of elevation and screening from 
nearby mountains at La Palma Island in Canary Island. Up- 
per left panel shows the terrain elevation in the island (lev- 
els indicate elevations of 0-500m, 500-1000m, 1000-1500m, 
1500-2000m, > 2000m), upper middle panel shows the com- 
posite radiance calibrated satellite image of the island, the 
upper right panel show the satellite image after deconvolu- 
tion. We superimposed the curves of equal elevation. Lower 
left panel shows the artificial sky brightness at sea level pre- 
dicted from the original satellite image. Lower middle panel 
shows the artificial sky brightness predicted from the de- 
convolved data when accounting for elevation. Lower right 
panel shows the predictions when accounting for elevation 
and mountain screening. The curves of equal elevation are 
superimposed. When accounting for elevation, the moun- 
tains at North and South of the sources diminish the sky 
brightness so that the isophotes appear more flattened. How- 
ever there are no effects where elevation is zero like e.g. in 
the sea near the top of the image. Screening by mountains 
has effects also at sea level and the 'umbra' of the moun- 
tains is visible looking carefully at the lower right panel. 
In particular, the screening by the Northern mountains is 



clearly visible near the top of the map. Colour levels from 
black to orange in the lower panels indicate zenith artificial 
brightness of <2.5, 2.5-5.0, 5-10, 10-20, 20-40, 40-80, 80-160, 
160-320, 320-640, >640 /icd/m 2 . Pollution at zenith is very 
low at the observatory site Roque de los Muchachos, well 
under the 10% of the natural sky brightness over which the 
International Astronomical Union consider the sky polluted. 
Maps refers to aerosol content K=l and do not accounts 
for local atmospheric conditions and denser aerosol layers. 
Results for Canary Island and Chile were not corrected to 
1998-1999 and refers to 1996-1997. 

Fig. |^ shows the effect of screening on the sky bright- 
ness in the nearby of Cerro Tololo Observatory for the same 
aerosol content. Colour levels indicate the same artificial 
brightnesses of fig. but we added two more levels. Left 
panel shows sea level sky brightness neglecting screening 
and a selected area. Right panel shows for the same area 
the contours when elevation and screening are accounted 
for, superimposed to the digital elevation map. La Serena is 
about 50 km far from Cerro Tololo so altitude and screening 
are much less effective than for La Palma. Their effects are 
recognizable only in the mountains surrounding the other 
two cities. 

Due to the larger computational time requested and to 
the fact that these effects are in general quite small in re- 
spect to the uncertainties, we neglected screening by terrain 
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elevation in the maps of section 
count whenever a greater precision is required and it become 
a normal practice in future when better code optimisation 
and faster workstations will be available. 



5 CONCLUSIONS 

We extended the method introduced in Paper I to map the 
naked eye star visibility and telescopic limiting magnitudes 
in large territories from DMSP satellite. This requires ac- 
counting for altitude of each land area, natural sky bright- 
ness in the chosen sky direction, stellar extinction in the cho- 
sen photometric band and eye capability. We also take into 
account mountain screening for near zenith sky directions. 
We presented, as an application, the maps of naked eye star 
visibility and total sky brightness in V band in Europe at 
zenith. Maps of limiting magnitudes in other directions will 
be useful to predict visibility of astronomical phenomena. A 
complete mappi ng of the brightness of the sky at a site, like 
Cinzano (20001), using satellite data instead of population 
data was already obtained (Cinzano & Elvidge, in prep.). 

Further improvements should be the development of: (i) 
a faster code to account for screening in a reasonable com- 
putational time, (ii) a set of worldwide three-dimensional 
atmospheric data sets or models for aerosol and molecules in 
order to change from standard clean atmosphere to the typ- 
ical night atmosphere in each territory in the given season, 

(iii) a method to measure the upward emission function of 
each land area from satellite data, solving problems like the 
presence of snow or fishing fleets (Cinzano et al., in prep.), 

(iv) large measurement campaigns in order to better con- 
strain the models (Cinzano & Falchi 2000). We hope that 
our work will be useful to the comprehension of how much 
mankind's perception of the Universe is endangered (see also 
Crawford 1991; Kovalevsky 1992; McNally 1994; Isobe & 
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Figure 7. The effects of elevation and screening from nearby mountains at La Palma Island. Upper panels show the terrain elevation 
(left), the composite satellite image (middle), the deconvolved satellite image (right). Lower panels show the artificial sky brightness at 
sea level (left) , accounting for elevation (middle), accounting for elevation and mountain screening (right). Superimposed are the curves 
of equal elevation. 




Figure 8. The effects of screening in the nearby of Cerro Tololo. Left panel shows sea level sky brightness. Right panel shows for the 
selected area the corresponding contours when elevation and screening are accounted. The digital elevation map is superimposed. 
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Hirayama 1998; Sullivan & Cohen 2000; for a large reference 
list see Cinzano 1994) and to support the battle against light 
pollution carried on worldwide by the International Dark- 
Sky Association (see the Web Site www.darksky.org) . 
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